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12.1 Introduction 



In 1931, Onsager 1 described the set of general linear equations, often called the thermodynamic 
equations of motion, which gave a macroscopic description of flows of particles or of energy as being 
driven by various thermodynamic forces. He observed that for most processes of interest, the flow is 
directly proportional to the corresponding force that drives it. Common examples cited include heat 
transfer where the vector flow of energy is proportional to the temperature gradient, or Ohm's law in 
which the electrical current density is proportional to the voltage gradient. Onsager suggested that each 
of the flows, J ; , in a given steady-state system is proportional to the sum of all forces, Xy, acting on the 
system. This can be written generally as 



J- = 5>«;X; 



(12.1) 



In a three-flow system, the flows would be written as 

Ji = L n Xi + L l2 X 2 + I13X3 J 2 = L 2l X l + i-22^2 + £23X3 
J3 = £31X1 + -L32X2 + -L33X3 



(12.2) 



These macroscopic quantities can be related to the microscopic properties of the materials through the L« 
parameters. In the case of thermoelectric (TE) materials, a two-flow system including electrical current 
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flow J, and heat flow q (as fluxes), and the forces of potential gradient W, and temperature gradient VT, 

or 

J = I U VV + L 12 VT q = L 21 VV + L 22 VT (12.3) 

Ohm's law is given by J = crE = — crW for an isothermal sample (VT — 0) with electrical conductivity a 
and the heat flux is q = —kVT for a sample with no electric fields present. For inhomogeneous samples 
or when temperature gradients are present, then variations in the chemical potential p, as well as the 
electrostatic potential V, must be included in the potential gradient. The forces and flows given in 
Equation 12.3 can be rewritten in terms of particle current density J. and entropy current density J 5 , as 

] q = -L qq VpZ - L qs VT J s = -L qs Wfl - L 5S VT (12.4) 

where p = p + eV is the electrochemical potential. The heat current density (or heat flux) and energy 
current density (or energy flux) are given by 

q=H s W=q+/IJ, (12.5) 

For charge carriers being electrons, e = —1.602 X 10~ 19 Coulombs, and 



(7 1 

/ = = T 

<79 2 2 

" e e p 



( L m L ss L qs \ K hL = S * = ea (12.6) 



where p is the electrical resistivity, k is the thermal conductivity, S " is the transport entropy per particle, 
and a is the TE power (or absolute Seebeck coefficient). Combining Equation 12.4, Equation 12.5, and 
Equation 12.6 gives 

J S = S%-|VT (12.7) 

W = (TS* + il)] q - kVT (12.8) 

V/Z= ~e 2 p] q - S*VT (12.9) 

Equation 12.8 shows the total energy current density has contributions from entropy and particle 
transport, while Equation 12.9 can be viewed as a generalized Ohm's law, which under isothermal 
(VT = 0) and chemically homogenous (V/I = eW) conditions reduces to the normal Ohm's law 
(VV=-ep] q ). 

There can be no accumulation of energy within an infinitesimal volume; therefore, the total energy 
current density W, given in Equation 12.8, must have zero divergence. For the one-dimensional case, this 
leads to Domenicali's equation for energy balance 

9 / 9T(x)\ 2 da(x) 

— Ux)— — ) = - P f+JT(x)— ^ (12.10) 

dx\ 6X ) 6X 

and the equation of heat flux 

q(x) = JT(x)a(x) - k(x)^- (12.11) 

dx 
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12.2 Modeling Techniques 



12.2.1 Exact 

An analytically exact technique has been described by Sherman, Heikes, and Ure~ for a p-n couple in the 
configuration described in Figure 12.1. 

In this technique the efficiency is determined by setting up the governing equations for the full 
generator, such that Equation 12.11 is written to include the n- and p-contributions as 



Q h = /r h [a p (r h ) - a n (r h )] - * p (r h ) 



dTp 
dx 



dr n 

-A n K n (T h )— — 
x=o dx 



(12.12) 



where the cross-sectional area of the p-leg is assumed equal to 1 cm and the cross-sectional area of the n 
leg is A n . In the actual device, the cross-sectional are of the p-leg is not, in general, equal to 1 cm ; 
however, if the ratio of cross-sectional areas between the p- and n-legs remains constant, then there is no 
change in the efficiency. Equation 12.12 can be inserted into the formula for thermodynamic efficiency 
17, = P /Qh, along with the output power, P = I 2 Ri, where the electrical current is given by the open 
circuit voltage divided by the total resistance, or 



(12.13) 



J r [a v (T) - a n (T)]dT 

The load resistance, R L , can be eliminated from the equation by rearranging Equation 12.13 as 

i\Rl + J o [p P CD + ^- IdxJ = f " [ap(D - On(T)]dT (12.14) 

and subtracting the necessary components to obtain the electrical output power as 

I 2 Rl = J f " [«p(T) - a n (r)]dr] - 7 2 [J o U(T) + ^]d*] 
The efficiency can then be written as 

^ [Op(T) - a n (T)]dT - /[ J o [p p (T) + ^- ldxl 



T h [ar,(T h ) - a„(r h )] 



k JT h ) dT 



dx | x =o 



A n K n (r h ) dr n I 

I dx \x=o 



(12.15) 



(12.16) 




-^x 



FIGURE 12.1 Single couple TE generator and coordinate system for the exact modeling technique. (Source: 
Sherman, B., Heikes, R.R., and Ure, R.W. Jr., /. Appl. Phys., 31(1), 1-16, 1960. With permission.) 
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Using, / p = I, and J n A n = -1 gives 



17, 



J^ h [op(r) - « n (r)]dr - ; p J q Pp (r)dx + ; n J q Pn (r>dx 



Th[ap(r h ) - a n (Th)] 



K p (T h ) dT 



/» 



dx 



K n {T h ) dT n 
x =o J n dx 



(12.17) 



This can be simplified through a change of variable: 

7n(T) = - — - r— and 

7n dT 



such that 



dy n (T) 
dT 



-r n (T) 



y v (T) 



y n (T) 



KpCO 

d^ 

Jp dr 



(12.18) 



(12.19) 



where r n (T) = T — ^ is the Thomson coefficient for the n-leg, and the electrical current density can 

, . dT 

be written as 



f h kJT) , f r >. cbc n 

J r c y n (T) J r c n dT 



dT 



(12.20) 



where L is the full length of the n leg. Using a similar procedure for the p-leg, the efficiency for the p-n 
couple is 



T?t 



Qh-Qc 



y P (r c ) - y n (r c ) + rja p (r c ) - « n (r c )] 
/ P <T h ) - /n<Th) + r h [a p (r h ) - « n (r h )] 



(12.21) 



In determining the maximum efficiency with this method, trial values for y„(T h ) and y n (T^) are 
approximated. Then/ p (T) and y n (T) are determined by numerical integration of Equation 12.18, and the 
efficiency 17, is calculated by Equation 12.21. This procedure is repeated for more trial values of y p (T h ) 
and y n (T h ) until a maximum in the efficiency is found. The values of L and J n are then determined by the 
functions y p (T) and y n (T) that yield the maximum efficiency, and the integration of Equation 12.20. 

The thermal efficiency given in Equation 12.17 for the module can also be written in terms of the 
efficiencies of the individual legs as 



a p (r)dr - ;, 



i?t 



L Pp(r)cb 



T c 



a n (r)dr + /, 



n Pn 

Jo 



(T)dx 



T h [aJT h )] 



K p (T h ) dT 



h 



dx 



T h [a n (T h )] 



K n (T h ) dT n 
h dx 



itHt) 



(12.22) 



where q h and q h are the heat flow densities (in W/cm ) into the hot-sides of the p-type and n-type legs, 
respectively. For the x direction defined in Figure 12.1, L = J, and J n A n = —I; therefore 

VQh p )+17n(Qh„) 



1?t 



Qk + Qh„ 



(12.23) 
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where Q h and Q h are the heat flows (in watts) into the hot-sides of the p-type and n-type legs, 
respectively. 

12.2.2 Iterative 

In the iterative technique ' a single leg (n-type or p-type) of the TE device is mathematically subdivided 
into many small segments such that there is a correspondingly small temperature gradient across each 
segment. In this condition, the material properties of a, p, and k for a given segment are assumed to be 
constant over the small temperature gradient across that segment. Equation 12.10 and Equation 12.11 
can be rearranged into a coupled set of first-order equations as 

dT(x) 1 



dx k(x) 

and 



[JT(x)a(x) - q(x)] (12.24) 



dfl(x) j Ja(x)q(x) 

-5p = p(x)/ 2 [l - Z(x)T(x)] - H ; (12.25) 

dx k(x) 

For dx— ► Equation 12.24 and Equation 12.25 can be written for the first segment as: 

T 1 = T +— UT a -q ] (12.26) 



1\ = la 



L/ 2 (l-^)-^^ldx (12.27) 



The output from one segment of the leg becomes the input for the adjacent segment, such that an 
iterative method can be used for determining the correct heat flow and temperature profile such as: 

1. Begin with a guess of q for the cold-side of the leg (at temperature T ). 

2. Use the values of X and q to calculate T Y and qi from Equation 12.26 and Equation 12.27, and 
the material properties of the leg at the temperature T . 

3. Repeat this procedure for calculating T 2 and q 2 from T Y and q u and the material properties of 
the leg at the temperature T ; . Continue for each segment up to the final segment at T N and q N . 

4. If T N is not equal to the hot-side temperature, then chose a new guess value for q and start the 
procedure over again. 

5. Continue this iterative process until T N is equal to the desired hot- side temperature. 

This is repeated for the other leg of the module. The efficiency of the n-type leg can then be 
determined as 



«nW — T dx + / n p(x)dx J 



fL 
In 

r, n = _«»_ -^ -^- -!- (12.28) 

IK 



and similarly for the p-type leg (under normal operation q a , q hn , and L are negative quantities for 
the x direction defined in Figure 12.2). 

To determine the maximum module efficiency based on r) n and i7 p , Equation 12.23 is maximized. 
This iterative technique was tested for the examples given in Ref. [3], with material properties 
described in Table 12.1. The module for Example I is for a temperature gradient of T c =400 K and 
T h = 750 K, while Examples II and III (fictitious materials) have T c = 400 K and T b = 1500 K. 

The cross-sectional area of the p-type leg is assumed equal to 1 cm . The optimal cross-sectional 
area for the n-type leg is then determined as described below, such that the ratio of areas is known. 
The cross-sectional areas can be varied without degradation to the maximum efficiency as long as this 
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FIGURE 12.2 Mathematical segmentation of each leg of a TE module. Equation 12.26 and Equation 12.27 are 
applied to each segment and iteratively solved from the cold to the hot end until the desired temperature profile is 
obtained and the heat flow determined. 



TABLE 12.1 TE Properties for Example Modules Described in Ref. [3] 



Example I (n-type) 
Example I (p-type) 
Example II (n-type) 
Example II (p-type) 
Example HI (n-type) 
Example III (p-type) 



<t (S/cm) 




a (/iV/K) 




k (W/cm K) 


T- 310 

O n = 

0.1746 


"n 


= 0.268-T- 


- 329 


K„ = 


54 
T 
3.194 

7^" 

T 
10 

Y 


o- p = 25 
10 5 

cr p = T 




= 0.150-7 + 211 
= 0.20-T - 400 


K P = 
K„ = 


« P 


= 200 




«p = 


o„ = 1000 


«„ 


= 0.20-T 




K n = 


3 
T 


o-p = T 


«p 


= 200 




K P = 


10 

Y 



Current Density for n-type leg (A/cm 2 ) 
20 25 30 35 40 45 50 



3 - 



2.5 - 



2 - 



1.5 - 



1 - 



£ 0.5 



- 



-0.5 




2.5 



j — i — I — i — i — i — i — I — i — i — i — i — I — i — i — i — i — I — i — i — i — i — I — i — i — i — «—J 2 4 



2.9 



2.8 



2-7 " 



2.6 



3= 
LU 



-2.5 -2 -1.5 -1 -0.5 

Current Density for p-type leg (A/cm 2 ) 



0.5 



FIGURE 12.3 Efficiencies for the n-type and p-type legs of Example I. 
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2.95 



2.9 



E 2 - 85 

o 

n 
o 



2.8 



2.75 




-1.6 



-1.4 -1.2 

^ (A/cm 2 ) 



-« — A„ = 0.029 (cm) 
■* — A n = 0.031 (cm) 
S—A„ = 0.032 (cm) 
/I n = 0.034 (cm) 



-♦— /I n = 0.035 (cm) 
-• — A„ = 0.036 (cm) 
-A- A„ = 0.038 (cm) 
-A- /I n = 0.039 (cm) 
-X- /l n = 0.041 (cm) 
-* — A n = 0.042 (cm) 
-+— A„ = 0.044 (cm) 



-I i i i_ 



-0.8 



FIGURE 12.4 Efficiency vs. L for various A n values for Example I module. 



ratio remains constant. Results for Example I materials give individual leg efficiencies as shown in Figure 
12.3, and Figure 12.4. For Example II there is a larger shift from the peak of rj p for the maximum 
efficiency of the module r/,, as shown in Figure 12.5, Figure 12.6, and Figure 12.7. The temperature 
profiles along the leg and heat flux through the n-type and p-type legs obtained from the iterative 
technique is shown in Figure 12.8. 



Current Density for n-type leg (A/cm 2 ) 
15 20 




27.8 



- 22 



_l I I I I I l_ 



14 



I 12 



24 



20 cd 



CD 

16 i 



-72 -68 -64 -60 -56 -52 -48 -44 -40 
Current Density for p-type leg (A/cm 2 ) 



FIGURE 12.5 Efficiencies for the n-type and p-type legs of Example II. 
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FIGURE 12.6 Efficiency vs. L for various A n values for Example II module. 

A comparison between results from the iterative technique and the exact approach for the three 
examples given in Ref. [3] are listed in Table 12.2. 

Agreement between the modeling techniques is found to improve with finer current density steps and 
smaller values of dx in Equation 12.26 and Equation 12.27. 

The iterative technique is also applicable to modules with segmented legs consisting of multiple n-type 
materials and/or multiple p-type materials. A test example found in Ref. [7] was used for a segmented 
structure. In this reference, the module efficiency is found by using the average material properties over 
the temperature gradient of a specific segment of the leg. The module consists of a p-type leg with three 
segments (P 1 , P 2 , and P 3 ) and an n-type leg with two segments (N l and N 2 ). The material properties are 
taken as temperature independent with values given in Table 12.3. 



25.85 




25.55 



5 5.5 
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FIGURE 12.7 Module thermal efficiency, 17,, vs. A n for Example II module. 



) 2006 by Taylor & Francis Group, LLC 



Modeling and Characterization of Power Generation Modules Based on Bulk Materials 12-9 



CD 



CD 
Q_ 

E 

CD 



(a) 



-20 



-22 



o) -24 

CD 



CO 
CD 



-26 - 



-28 



0.2 0.4 0.6 0.8 

Position along the leg (cm) 



1500 




i — i 


n-rrn 




1 — i — I 








n 


r" ~ 














-type 






1400 


F 






— • — n-type 
































i_ 






- 


1300 

1200 

1100 

1000 

900 

800 

700 

600 

500 

400 




i_l 


l_i 




1 


l_i_i 








5.5 



- -6 



en 

CD 



-6.5 £ 



(b) 



0.2 0.4 0.6 0.8 1 

Position along the leg (cm) 



FIGURE 12.8 Temperature profiles (a) and heat flux (b), q n and q p , on the n-type and p-type legs for the Example II 
module. 



A comparison of the results from the iterative technique and the results given in Ref. [7] are shown in 
Table 12.4. 

The dimensions of each segment are determined from plots of the material properties for each leg such 
as shown in Figure 12.9. 

12.2.3 Averaging Technique 

A simple approach for estimating the efficiency can be obtained by using average values for the material 
properties of the module legs. In this situation, the Thomson coefficient is equal to zero, and the 
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TABLE 12.2 Comparison between Iterative and Exact Modeling Techniques 





Iterative Technique 


Exact Approach 


Example I 






Efficiency (%) 


2.94 


2.9 


Current (A) 


1.24 


1.22 


A n (cm 2 ) 


0.0365 


0.036 


Example II 






Efficiency (%) 


25.8 


26 


Current (A) 


55 


55 


A n (cm 2 ) 


4.48 


4.5 


Example III 






Efficiency (%) 


33.5 


31 


Current (A) 


65.2 


65 


A n (cm 2 ) 


2.20 


2.2 



TABLE 12.3 Segmented Module from Ref. [7] 



Segment 


u (S/cm) 


k (W/m K) 


a (rxV/K) 


Temperature Span (K) 


Pi 


1886.8 


4.636 


154 


720 -880 


Pi 


277.0 


1.258 


226 


420 -720 


P 3 


1357.4 


1.861 


159 


300 -420 


Ni 


246.9 


1.342 


231 


420 -880 


N 2 


769.23 


1.453 


212 


300 -420 



efficiency is given by (for A p = 1 cm ) 



Vi 



I 2 R 



n )T h I+ ^(k v +A n K n ) - I 2 l\(p? + f 1 ) 



O p - a 



(12.29) 



TE devices are heat pumps and are Carnot efficiency-limited. Equation 12.29 can be rearranged as the 
product of the Carnot efficiency and the thermocouple efficiency as 



li 



T h - t c Vi + zr - i 



h Vi+zf- 



r„ 



(12.30) 



TABLE 12.4 Results from the Test Structure Described in Ref. [7] 



Length of Pi (mm) 
Length of P 2 (mm) 
Length of P 3 (mm) 
Length of N; (mm) 
Length of N 2 (mm) 
A n (cm 2 ) 

■n (%) 

1(A) 



From Ref. [7] 



Iterative Technique 



5.51 

2.91 

1.58 

7.72 

2.28 

1.985 

11.48 

32.135 



5.46 

2.99 

1.57 

7.93 

2.07 

1.964 

11.48 

32.60 
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FIGURE 12.9 Profile of the (a) thermopower for the segmented p-type leg and (h) electrical conductivity for the 
n-type leg found from the iterative technique. 

1 



where T = — (T h + T c ) and 



(" P ~ «n) 
L Vp K P + %/PnKn J 



(12.31) 



is the figure-of-merit for the thermocouple. The efficiency is maximized for 



,„=l^ and* = 7171^+^ 



or the load resistance is equal to VI + TZ times the resistance of the module. 



(12.32) 
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12.3 Multidimensional Analysis 

The methods described in the previous sections assume the temperature and the electric potential 
distributions in each p- and n-leg of the TE module to be one-dimensional. Although this is generally a 
good assumption, there are configurations for which this assumption may not be valid. In this section, 
methods are presented that do not make this assumption and enable the analysis of legs, where the 
temperature and the electric potential distribution can be one-, two-, or three-dimensional. The 
organization of this section is as follows. First, the governing equations are described. Next, an efficient 
finite-volume (FV) method of solution is presented. This section concludes with computed results 
illustrating the multidimensional effects that can arise from heat transfer. 

12.3.1 Governing Equations 

The equations governing the multidimensional temperature and electrical potential distributions in TE 
materials under steady-state conditions and in the absence of an applied magnetic field are the energy 
equation and current flow equation given by Refs. [9-12]: 

V(Kvr) + P j 2 -rj|(^)vr + (Va) r ] = o (12.33) 



V-J = (12.34) 



where 



J= -0- V(- + Vj + aVr (12.35) 

q=aTJ-/<VT (12.36) 

As in previous sections, the vector J is the electric current per unit area; the vector q is the heat 
transfer rate per unit area; T is temperature; k is thermal conductivity at zero current; a = 1/p is 
electrical conductivity; p is electrical resistivity; a is the absolute Seebeck coefficient (thermopower); 
fi is the chemical potential; V is the electrostatic potential; and e is the charge of charged particles 
giving rise to the electrical current. Note that a, k, a, and p are functions of temperature. 

Equation 12.33 to Equation 12.36 form a system of two coupled partial differential equations (PDEs) 
with two dependent variables, the temperature T and the electrical potential V. The boundary 
conditions (BCs) for the temperature are as follows. At all surfaces exposed to vacuum or gas, the BC 
imposed is either specified heat flux or specified heat transfer coefficient, i.e. 

q" = q's (12.37a) 



1 1 convection ' Eradiation (11.5/0) 

where 

Hconvection "\ -* s -* 00/ \\-Z.j / C) 



II 



q'radiation = ecr'T* + (1 - a op )G s ~ h md (T s - T x ) (12.37d) 

In Equation 12.37a to d, T s and q" are, respectively, the temperature and heat flux at the surface of 
the leg exposed to vacuum or gas; T s and q" can be constant or function of position along the 
surface; T^ is the bulk mean temperature of the gas in the cavity between the p- and n-leg; h is the 
heat transfer coefficient and its value depends on the mode of convection such as natural or forced; 
e and a op are, respectively, the emissivity and the optical absorptivity of the leg surface; a ' is the 
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Stefan -Boltzmann constant; and h mA is an effective radiation heat transfer coefficient. As written, 
Equation 12.37d neglects gas radiation, which should be modeled if there is water vapor or carbon 
dioxide in the gas. At all surfaces of the leg that are connected to the conducting plates, the BC 
imposed is either specified temperature or heat flux, i.e. 

T = T Q (12.38a) 



q = q (12.38b) 

where T and q can be constant or function of position along the surface. 

The BCs for the electrical potential are as follows. At all surfaces of the leg that are exposed to vacuum 
or gas, the current must be parallel to the surface, i.e. 

J-n=0 (12.39) 

where n is a unit vector normal to the surface. At all surfaces of the leg that are connected to 
the conducting plates, the BC imposed is either specified electric potential or continuity of electric 
current, i.e. 

V = V (12.40a) 



J = Jo (12.40b) 

In Equation 12.40a and b, V a and J can be constant or function of position along the surface. 

12.3.2 Numerical Method of Solution 

Since the PDEs governing the temperature and electric potential distributions in TE materials given by 
Equation 12.33 to Equation 12.36 are not linear and the BCs given by Equation 12.37a to d to Equation 
12.40a and b can be complicated functions of temperature and position, it is difficult to obtain exact 
solutions. However, there are a wide variety of numerical methods that can be used to generate 
approximate solutions to these equations (see for example, Refs. [13,14]). In this section, an efficient FV 
method that can handle complicated geometric shapes for TE legs and modules is presented. 

Before presenting the FV method, it is noted that when solving Equation 12.33 to Equation 12.40a and b 
by any FV method, the FV method can be applied directly to those equations, or to an unsteady or time- 
dependent form of those equations. If a FV method is applied directly to Equation 12.33 to Equation 12.40a 
and b, then a system of nonlinear algebraic equations is obtained. Though this system can be solved 
iteratively by using the Newton -Raphson method or one of its variations, convergence requires the initial 
guess to be within the radius of convergence, " which is a difficult task. By applying a FV method to an 
unsteady form of the governing equations, initial guess can be constructed from physical considerations, 
and convergence is achieved when the time-derivative terms added to the governing equations approach 
zero within some prescribed tolerance. In addition to an easier path to convergence, the latter approach 
offers greater flexibility in constructing efficient algorithms. Thus, Equation 12.33 and Equation 12.34 are 
cast in unsteady form as follows: 



i ar 

y r dt 



V(kVT) +/(r, n /(r,io = pj 2 -rj|7^)vr + (va) T l (12.41) 

1 dV 
— — = -V-J (12.42) 

y v dt 

In the above equation, t denotes time or pseudo-time. Since time-accurate solutions are not of 
interest, y T and y v can take on any value that accelerates convergence to steady-state. By adding 
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a time-derivative term to Equation 12.33 and Equation 12.34, these equations are converted from 
elliptic to parabolic, but they revert back to elliptic once the time -derivative terms on the left-hand 
side of Equation 12.41 and Equation 12.42 vanish. Since the coefficient of the second derivative 
with respect to space coordinates must be positive for well-posed parabolic PDEs, the sign of y T 
and y v must be positive. 

There are many FV methods that can be used to generate steady-state solutions to Equation 12.41 and 
Equation 12.42. These include explicit, implicit, and hybrid methods. Explicit methods (e.g., methods 
in which the time derivatives are approximated by Runge-Kutta or Adams -Bashforth methods ) are 
extremely easy to derive and program, but they are conditionally stable. For the energy equation given by 
Equation 12.41 explicit methods can be used, because the stability criteria do not impose overly small 
time-step sizes so that convergence to the steady-state solution can be obtained with reasonable 
efficiency. However, for the electric-current equation given by Equation 12.42, the use of explicit methods 
does impose a very small time-step size, which will greatly slow convergence to steady-state. Since 
Equation 12.41 and Equation 12.42 are coupled, the smaller time-step size imposed by Equation 12.42 
must be used. Thus, explicit methods are inefficient for solving Equation 12.41 and Equation 12.42. 
Implicit methods are much more stable than explicit methods; they are unconditionally stable for linear 
PDEs. By allowing for orders of magnitude larger time-step sizes, implicit methods converge much faster 
than explicit methods. However, implicit methods are more complicated to derive and program. In this 
study, a hybrid method is selected in which Equation 12.41 and Equation 12.42 are solved by using an 
explicit and implicit method, respectively. With the hybrid method, the maximum time-step that can be 
used is limited by the explicit method applied to Equation 12.41. 

To illustrate the hybrid method for obtaining solutions to Equation 12.41 and Equation 12.42, consider 
a leg of a TE module with length L x , and a rectangular cross-section ofL, byi z (Figure 12.10). As with all 
FV methods, the derivation of the hybrid method for Equation 12.41 and Equation 12.42 along with BCs 
given by Equation 12.37a to d to Equation 12.40a and b involve the following three major steps: 

(1) Discretize the domain; 

(2) Replace the PDEs by algebraic equations; 

(3) Specify the solution algorithm. 

Each of these three steps is described below. 

Discretize Domain: In this study, the domain of the problem (the rectangular leg shown in Figure 12.10) 
is discretized by dividing it into (IL — 1)(JL — 1)(KL — 1) identical size cells as shown in Figure 12.11. 
The center of each of these cells is located at 



(i-i)Ax, Ax = L x /(iL - 1), i = 1,2, ...,IL — 1 



(12.43a) 




FIGURE 12.10 Leg of a TE module. 
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FIGURE 12.11 Leg divided into cells. 



yj = (j-l)Ay, Ay = L/(JL - 1), j= 1,2 JL- 1 (12.43b) 

Zk = (k - I)Az, Az = L Z /(KL - 1), fc= 1,2,...,KL- 1 (12.43c) 

and each cell is bounded by the following coordinates: 

Xi-m and x i+m , y hm and y ; - +1/2 , z t _ 1/2 and z k+1/2 (12.44) 

In the above, IL, JL, and KL are integers, and their values determine the number and the shape of the cells 
in the domain. With such a simple discretization, the volume of each cell is 



vyjt = AxAyAz 
and the areas of the faces bounding each cell are 



A i~V2,j,k 



HjM-m 



H+\n,j,k 



Hj+l/2,/c 



A 



i,j,k+l/2 



AyAz 
AxAz 
AxAy 



(12.45a) 

(12.45b) 
(12.45c) 
(12.45d) 



If the cells in the domain are of unequal sizes, then Equation 12.45a to d must be modified to reflect the 
actual volumes and surface areas. 

Discretize PDEs: With the domain discretized, the next step is to replace the governing equations by 
algebraic equations. This is accomplished in three steps. The first is to integrate the governing equations 
about each cell. For the (i,j, fc)fh cell, integration of the energy and electric-current equations given by 
Equation 12.41 and Equation 12.42 yields: 



\Av T ^-^ T) - f(T > v) ) dv 

i.j.k 

[ {— ^+V-j)dv = 
Jv, M . V 7v 9f / 







Upon noting that none of the cells deforms in time, the above equations become 



— (tt) =^)«.M+f V(KVT)dv 
J T \ 9f /i,j,k Jv,, M . 



(12.46) 
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3*(£) =[ "V.Jdv (12.47) 

Tij, k = — f rdv, Vy,t=— f vdv, /y*= — f /(r,v)dv (12 



48) 

where T and V are, respectively, the average temperature and electric potential in the (/', j, fc)th cell. By 
invoking Gauss' theorem, the volume integrals involving divergence can be converted to surface integrals 
so that Equation 12.46 and Equation 12.47 become 



v 
Jt 



— (lP\ =Cfv) i , j , k +[ (KVT)-ndA (12.49) 

Yt \ 9^ /ij,k J„„ 

^(£) = [ J-ndA (12.50) 

n = n x i + nj + n z k (12.51) 

In the above equations, 9v denotes the surface bounding the cell; n denotes the unit vector normal to the 
cell surface and pointing away from the cell; and n x , n ; „ and n z are directional cosines. Since each cell in 
Figure 12.11 are bounded by six faces, and the cell-face normals are aligned with the x, y, or z coordinate 
direction, Equation 12.49 and Equation 12.50 become 

V / rt T \ /face6 \ 

fUlr^Hl^H,, <1252a) 

Jv V 81 ):# Vtoi Aji 



where 

facc6 



iacc6 

X ^' n &ce =(/* A )j+ijJt ~ (/x^-);-l,j,t + (7 r A)y+i,/t ~ (VOy-l./t 
facel (12.53b) 

+ (JzA) i#+ l -(/z^)y ; it-I 

; *=-'{^(7 + v ) + a ^} ^ = ^< z (12 ' 53c) 

In the above equations, the bars over J x ,J y ,J z , and the derivatives of T denote averages over the respective 
cell faces. In this study, all derivatives in the fluxes across the faces as well as those appearing in the source 
term /are approximated by second-order accurate central differences. The fluxes in the x direction are 
approximated as follows 






(*A),. +Ut T ' +m . T '' hk (12.52c) 
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(*A)«-ij,* 



CIA) 



x A )i+lj,k 



"(«A) 



(£ + ?) -U + r) 

\ e )i+i,j,k \ e /,j,k 



Ti,j,k Ti-l,j,k 



= -{#£), 



Ax 



Ax 



(12.52d) 



a 



i+ij,k- 



,i,k\ 8 xi — + v ) L . + a i+{,j,A T i+{,j,k 



Ti+l,j,k Ti,j,k 

Ax 



(12.53d) 



ClA\-l, h k= -(^)i-i,;,/ 



\ e )i,j,k \ e /i-i,j,k 



Ax 



"■{- 1 ik~ 

1 2 V' K 



Ax 



-(aA)i- 



4 s ^ +v U, 



- a i-±,j,k S xTi-} 



,i,k 



(12.53e) 



$i±Uk = l(0i±l# + <A;,;,/tX i[i= K,a,a 



(12.54) 



The derivatives of V, T, and /u, in the source terms are approximated as follows: 



V ox Ji,j,k 



$, 



i+Uk 



(p 



i-hJ- k 



Ax 



<P 



I— 2 J,/C 



®i±Uk + ®i,j,h <P=V,T,il 



(12.55) 



The fluxes and derivatives in the y direction and z direction are approximated in a similar manner. 

For all cell faces that are on the boundary of the domain, the appropriate BC given by Equation 12.37a 
to d to Equation 12.40a and b must be applied. For example, the heat flux at all surfaces of the TE leg 
exposed to gas must satisfy Equation 12.37a to d. Also, the electric potential there must be such that the 
electric current only flows along the surface with no component perpendicular to the surface as 
stipulated by Equation 12.39. 

The time derivative in Equation 12.52a is approximated by the conditionally stable, first-order 
accurate in time, Euler explicit method. Since only the final steady-state solution is of interest, the 
lowest-order accurate time-difference formula is used. With the Euler explicit method, Equation 12.52a 
becomes 



1 i,j,k — i i,j,k " 



y T At(fl k ) + ^ ( f 6 *Vr-n face ) 



(12.56) 



In the above equation, the superscript n denotes data at the nth time level (i.e., t" = nAt; At is taken to 
be a constant), where the solution is assumed to be known, and the superscript (n + 1) denotes data at 
the (n + l)th time level (i.e., t" +1 = t" + At), where the solution is unknown and to be determined. 
Equation 12.56 is the algebraic representation of the energy equation given by Equation 12.41, and it is 
valid for all cells whose faces do not coincide with a boundary. For those cells with faces that coincide 
with a boundary, appropriate BC for the temperature or heat flux must be imposed. Thus, we rewrite 
Equation 12.56 as follows so that it is valid for all cells: 



ffn+l 



1 i,j,k ' 



„ », /face6 \ " 

V, j> k Vfacel /,',; 



(12.57a) 
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where 



q = — kVT for all faces that do not see gas (12.57b) 



1 = Iconvection + Eradiation P er Equation 12.37a to d for faces that see gas (12.57c) 

For Equation 12.53a the first-order accurate in time, Euler implicit method is used to approximate the 
time derivative term, which is numerically more stable than the Euler explicit method. Again, only the 
lowest-order accurate in time method is used because only the steady-state solution is of interest. By 
using the Euler implicit method, Equation 12.53a becomes 

^-nUt-V^IIJ-"*-) (12.58) 

»j.* Vfacel / i,j,k 

Unlike Equation 12.56, the right-hand side of Equation 12.58 involves derivatives of T and V at the 
(n + l)th time level, where they are unknown instead of the nth time level, where they are known. The 
unknown f" +1 in Equation 12.58 can be obtained from Equation 12.57a to c to be explained later, but 
the unknown V" +1 requires the solution of a system on linear equations of order m = (IL — 1) 
(JL - 1)(KL - 1) and of bandwidth w= min[(IL - 1)(JL - 1), (IL - 1)(KL - 1), (JL - 1)(KL - 1)]. 
The operation count needed to solve this system of equations can be reduced by orders of magnitude if a 
method known as approximate factorization (AF) is used. 14 To implement the AF method, Equation 12.58 
is rewritten in the following operator form: 



where 



0AVl k + (<fe + 4> y + <t> z )AVlj,k = 4>r (12.59a) 

0= _h>±_ (12.59b) 
y v At 

^li,k = V"u ~ Vlj,k (12.59c) 

<t> x = -[( < 7AS x ) !+ i^-( ( fAS ;c )i_i # ] (12.59d) 

^ = -[fa^A^ 8 ?)^] (12 - 59e) 

4 = -[(aA8 z ) tj ^-(aA8 z ) ljk _ { ] (12.59f) 

$t = ~(A + <t> y + <t>J ~ + V*J -(4>'x + <& + &)Tlt k l (12.59g) 

4>' x =-[(aaAd x ) i+hJtk -(ffaA8 x )^ Uk ] (12.59h) 

4>' y = -\{ffaA8 y )., +kk -{ffaA8 y ).._ k ^ (12.591) 

<t>' z = -[(ffaA8 z ) ijtk+i -{aaA8 z \ hk _ , ] (12.59J) 
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Equation 12.59a can be approximately factored as follows 17 : 

(6+ <i> x )B-\6+ <t> y )6-\e+ jJAVfa = 4> T (12.60) 

Note that at steady-state, AV" k = so that the errors introduced by the AF vanish. Equation 12.60 can be 
split into the following three one-dimensional operators: 

{0+<f> x )AVl JM =<i>T (12.61a) 

(6+cl> y )AV*l k =0AVl Jtk (12.61b) 

{6+<l> z )AVl hk =0AV** x (12.61c) 

V"£ = V" m + AVl k (12.61d) 

Equation 12.61a to d forms systems of linear equations with tridiagonal coefficient matrices that can be 
solved very efficiently by the Thomas algorithm. 14 Also, the order of the matrices is now (IL — 1),(JL — 1), 
or (KL - 1) instead of [(IL - 1)(JL - 1)(KL - 1)]. 

Specify Solution Algorithm: With the domain discretized and the governing equations replaced by 
algebraic equations, the solution algorithm can now be described as follows: 

1. Input geometry of leg by specifying L x , L y , and L z . 

2. Input TE material properties by specifying k, a, a, and p as a function of temperature. 

3. Specify BCs for T and V via Equation 12.37a to d to Equation 12.40a and b. 

4. Calculate cell volumes, cell face areas, and spacing between cell centers by using Equation 12.43a 
to c and Equation 12.45a to d. 

5. Specify y T , y v , and the time-step size (Af) to be used that ensure numerical stability, and rapid 
convergence to steady-state solutions. 

6. Set the time level index n — 0. 

7. Specify the initial condition or initial guess for T and V at every cell (i.e., T", k and V"- k at 
n = 0). T" ■(. and V": k at n = can be constant or linearly interpolated from their values on the 
boundaries. 

8. Calculate TE material properties (k, a, a, and p) at every cell based on T": k . 

9. Calculate f at every cell at the (n + l)th time level {T"l~ k ) by using Equation 12.57a to c. 

10. Calculate TE material properties (k, a, a, and p) at every cell based on T"/fc ■ 

11. Calculate Vat every cell at the (n + l)th time level (V"~f~ k ) by using Equation 12.61a to d, where 
T?f. x is from Step nine. i~n+i _ r" I \V" +1 — V" I 

, ^ ,, • „ I 1 i.i.k 1 i,j,k\ . , \ i.j.k i,j.k\ n .1 .1 

12. Examine convergence. It r= — r~r, < s r and r-= — ~r, < s v at every cell, then the 

1% I \ V l,,k\ 

solution at the in + l)th time level is taken to be converged. Otherwise, set n to n + 1 and repeat 
Steps 9 to 12. In the above, the convergence criteria, s T and e v , can be set as 10 5 or some 
smaller number. 

12.3.3 Computed Results and Discussions 

In this section, the numerical method presented in the above is applied to study the temperature and 
electric potential distributions in the leg of a TE module depicted in Figure 12.12. Table 12.5 summarizes 
the four cases simulated to illustrate how heat transfer can produce multidimensional effects. As noted in 
Figure 12.12 and Table 12.5, the TE leg investigated is rectangular in shape and has length L x (10 mm) 
and cross-sectional area L y L z (25 mm ). The temperature at the hot end of the leg (x = 0) is maintained 
at T h — AT(y/L y ), where T h is 900 K and AT is either K or 50 K. AT is intended to examine the effects 
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0.15 



0.1 



0.05 



h = 

h = 100W/m 2 -K 

h = 200 W/m a -K 

h = 200 W/m a -K 
AT=50K 




0.008 0.01 



h = 
h = 100W/m 2 -K 
h = 200 W/m 2 -K 

h = 200 W/m a -K 



(b) 




AT=50K 



0.002 



0.004 0.006 
x(m) 



0.008 0.01 



FIGURE 12.12 (See color insert following page 31-12.) Predicted Tand V distributions along x in an x—y plane 
midway between the ends in z direction. Multiple lines indicate variations in Tand V in the x-y plane. 



of nonuniform temperature distribution at the base of the TE leg, which can be appreciable if the 
cross-sectional area of the leg is large and if the heat transfer rate to or from the cold and hot sources are 
appreciable. The temperature at the cold end of the TE leg (x = L x ) is maintained at T c , which was set to 
be a constant at the ambient temperature of 300 K. The temperature BC for all surfaces exposed to the gas 
was specified by using Equation 12.37a to d, where h is either W/m K, 100 W/m K, or 200 W/m~ K, 
representing contributions from both natural convection and radiation heat transfer. T^ is set to be the 
arithmetic average of T^ and T c . The BCs used for the electric potential are as follows. At the hot end of 
the TE leg (x = 0), current is extrapolated, whereas at the cold end (x = L x ), the electric potential is set to 
zero. On all surfaces of the TE leg exposed to the gas, Equation 12.39 is imposed. In all cases, the TE 



TABLE 12.5 Summary of Cases Simulated 



Case 


7h(K) 


r c (K) 


AT(K) 


h (W/m 2 


K) 


L x , L„, L z (mm) 


1 


900 


300 










10, 5, 5 


2 


900 


300 





100 




10, 5, 5 


3 


900 


300 





200 




10, 5, 5 


4 


900 


300 


50 


200 




10, 5, 5 
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Case 3 



Case 4 



900 

BOO 
700 
600 
500 
400 
300 




FIGURE 12.13 (See color insert following page 31-12.) Predicted Tand V distributions in the middle x—y plane. 

material (n-type) has the following temperature-dependent properties 

a = -297.7 + 1.855T - 0.0061604T 2 + 6.4013 X KT 6 T 3 - 2.1847 X KT 9 T 4 (fiV/K) 

a= 116.83 + 13561 e " 00068336T (S/cm) 
k = 0.23912 + 609.95/T - 838.3/T 2 (W/m K) 

The results generated for Cases 1 to 4 are plotted in Figure 12.12 through Figure 12.14. In Figure 12.12, it 
can be seen that for Case 1 with h = (i.e., no heat loss from the four sides of the TE leg), the 
temperature profile is nearly linear. This is because, the effect of having a thermal conductivity that 
decreases with temperature is offset by temperature increases from heat generation due to electric current 
flow. For Case 1, Tand Vare onedimensional (i.e., they are only functions of x), and serve as a reference 
for comparing with cases for which there are heat transfer from the four sides of the TE leg exposed to gas. 

For Case 2 with h = 100 W/m K, it can be seen from Figure 12.12 and Figure 12.13 that the Tand V 
become three-dimensional (3-D). 

Since T x = (T h + T c )/2, the TE material losses "heat" when its surface temperature is higher than T&, 
and gains "heat" when its surface temperature is lower than T m . Thus, temperature near the hot end is 
lower than the 1-D prediction, whereas the temperature near the cool end is higher than the 1-D 



Case 3 



Case 4 




j J(A/m d ) 

120000 

00000 
80000 
60000 



FIGURE 12.14 (See color insert following page 31-12.) Predicted magnitude of the electric current in the middle 
x—y plane (top) and in a plane near the surface exposed to gas (bottom). 
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prediction. This observation can be discerned in Figure 12.12 by comparing the results obtained for Cases 
1 and 2. Once the temperature becomes multidimensional, Figure 12.12 and Figure 12.13 show that the 
electric potential distribution becomes multidimensional. By increasing the heat transfer to- or from 
the four sides of the TE leg exposed to gas (e.g., from h = 100 W/m K of Case 2 to h = 200 W/m K of 
Case 3), Figure 12.12 shows the temperature near the hot end to be lowered more and the temperature 
near the cold end to be increased more. In Case 4, the hot end is assumed to have a linear temperature 
variation of 50 K, which produces a nonsymmetric distribution of T and V. The effect of the 
multidimensionality in Tand Vis to make the electric current flow to be nonuniform throughout the TE 
material as shown in Figure 12.14. Not shown here is that when the gradients of Tand V have opposite 
signs, the distribution of Tand V differ markedly from those shown in Figure 12.12. 



Nomenclature 

Symbol Quantity 



Symbol Quantity 



A n 


cross-sectional area of the n-type leg 


T h 


temperature on the hot-side of the 


A p 


cross-sectional area of the p-type leg 




module 


E 


electronic charge (—1.602 X 10 C) 


T n 


temperature of the n-type leg 


E 


electric field 


T 


temperature of the p-type leg 


h 


heat transfer rate 


V 


electrostatic potential 


1 


electrical current 


W 


energy flux 


J 


electrical current density 


x, 


force giving rise to the flows 


J, 


flow of particles or energy 


a 


thermoelectric (TE) power 


\ 


particle current density 


Vn 


efficiency of the n-type leg 


\ 


entropy current density 


Vp 


efficiency of the p-type leg 


h 


Onsager parameters linking forces and 


Vt 


thermodynamic efficiency 




flows 


K 


thermal conductivity 


Po 


output electrical power from the module 


K n 


thermal conductivity of the n-type leg 


q 


heat flux 


K P 


thermal conductivity of the p-type leg 


q" 


surface heat flux 


fl 


chemical potential 


Qc 


heat flow on the cold-side of the module 


M 


electrochemical potential 


Q h 


heat flow on the hot-side of the module 


P 


electrical resistivity 


Rl 


electrical load resistance 


l"n 


Thomson coefficient of the 


s ,:: 


transport entropy per particle 




n-type leg 


T 


temperature 


T P 


Thomson coefficient of the 


T c 


temperature on the cold-side of the 




p-type leg 




module 


x,y,z 


position variables 
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